Problema 1

La siguiente figura sugiere como estimar el valor de \(\pi\) con una simulación. En la figura, un circuito con un área igual a \(\pi/4\), está inscrito en un cuadrado cuya área es igual a 1. Se elige de forma aleatoria \(n\) puntos dentro del cuadrado. La probabilidad de que un punto esté dentro del círculo es igual a la fracción del área del cuadrado que abarca a este, la cual es \(\pi/4\). Por tanto, se puede estimar el valor de \(\pi/4\) al contar el número de puntos dentro del círculo, para obtener la estimación de \(\pi/4\). De este último resultado se encontrar una aproximación para el valor de \(\pi\).

Imagen con circulo circunscrito
Imagen con circulo circunscrito

Pasos sugeridos:

Genere n coordenadas \(x: X1, . . . , Xn\). Utilice la distribución uniforme con valor mínimo de 0 y valor máximo de 1. La distribución uniforme genera variables aleatorias que tienen la misma probabilidad de venir de cualquier parte del intervalo (0,1).

Genere 1000 coordenadas \(y: Y1,...,Yn\), utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.

set.seed(123)
n1 <- 1000
n2 <- 10000
n3 <- 100000
x_1000 <- runif(n1, 0, 1)
y_1000 <- runif(n1, 0, 1)

x_10000 <- runif(n2, 0, 1)
y_10000 <- runif(n2, 0, 1)

x_100000 <- runif(n3, 0, 1)
y_100000 <- runif(n3, 0, 1)

Cada punto (Xi,Yi) se encuentra dentro del círculo si su distancia desde el centro (0.5,0.5) es menor a 0.5. Para cada par (Xi,Yi) determinesi la distancia desde el centro es menor a 0.5. Esto último se puede realizar al calcular el valor (Xi−0.5)2+(Yi−0.5)2, que es el cuadrado de la distancia, y al determinar si es menor que 0.25.

¿Cuántos de los puntos están dentro del círculo?

puntos_dentro_1000 <- sum((x_1000 - 0.5)^2 + (y_1000 - 0.5)^2 <= 0.25)
puntos_dentro_10000 <- sum((x_10000 - 0.5)^2 + (y_10000 - 0.5)^2 <= 0.25)
puntos_dentro_100000 <- sum((x_100000 - 0.5)^2 + (y_100000 - 0.5)^2 <= 0.25)

El numero de puntos dentro del circulo es de: 800

¿Cuál es su estimación de \(\pi\)?

\(\pi = 4 \times \frac{\text{puntos dentro}}{n}\)

pi_estimado_1000 <- 4 * puntos_dentro_1000 / n1
pi_estimado_10000 <- 4 * puntos_dentro_10000 / n2
pi_estimado_100000 <- 4 * puntos_dentro_100000 / n3

La estimación de \(\large\pi\) con \(\large n=1000\) es de: 3.2

La estimación de \(\large\pi\) con \(\large n=10^{4}\) es de: 3.1528

La estimación de \(\large\pi\) con \(\large n=10^{5}\) es de: 3.144

¿Cuál es el error absoluto de su estimación?

\(\Large|\pi - \pi_{estimado}|\)

error_absoluto_1000 <- abs(pi - pi_estimado_1000)
error_absoluto_10000 <- abs(pi - pi_estimado_10000)
error_absoluto_100000 <- abs(pi - pi_estimado_100000)

Error absoluto para muestra de \(n=1000\) : \(0.0584073\)

Error absoluto para muestra de \(n=10^{5}\) : \(0.0112073\)

Error absoluto para muestra de \(n=10^{4}\) : \(0.0024073\)

¿Cuál es el error relativo de su estimación?

\(\Large\frac{|\pi - \pi_{estimado}|}{\pi} = \frac{|\pi - 3.2|}{\pi}\)

error_relativo_1000 <- abs(pi - pi_estimado_1000) / pi
error_relativo_10000 <- abs(pi - pi_estimado_10000) / pi
error_relativo_100000 <- abs(pi - pi_estimado_100000) / pi

Error relativo para muestra de \(n=1000\) : \(0.0185916\)

Error relativo para muestra de \(n=10^{5}\) : \(0.0035674\)

Error relativo para muestra de \(n=10^{4}\) : \(7.6628216\times 10^{-4}\)

Problema 2

La simulación ayuda a entender y validad las propiedades de los estimadores estadísticos como son. insesgadez, eficiencia y la consistencia principalmente. El siguiente problema permite evidenciar las principales características de un grupo de estimadores propuestos para la estimación de un parámetro asociado a un modelo de probabilidad.

Sean \(X1, X2, X3 \,y\, X4\), una muestra aleatoria de tamaño \(n=4\) cuya población la conforma una distribución exponencial con parámetro \(\theta\) desconocido. Determine las características de cada uno de los siguientes estimadores propuestos:

  1. \(\hat{\theta}_1 = \frac{X_1 + X_2}{6} + \frac{X_3 + X_4}{3}\)

  2. \(\hat{\theta}_2 = \frac{X_1 + 2X_2 + 3X_3 + 4X_4}{5}\)

  3. \(\hat{\theta}_3 = \frac{X_1 + X_2 + X_3 + X_4}{4}\)

  4. \(\hat{\theta}_4 = \frac{\min\{X_1, X_2, X_3, X_4\} + \max\{X_1, X_2, X_3, X_4\}}{2}\)

n <- 4
m = 5000
lam <- 0.5

rep = rexp(m*n,rate = lam)
df = data.frame(matrix(rep,nrow = m, ncol = n, byrow = TRUE))

estimador1 <- function(m) { ((m[1] + m[2]) / 6) + ((m[3] + m[4]) / 3) }
estimador2 <- function(m) { (m[1] + 2*m[2] + 3*m[3] + 4*m[4]) / 5 }
estimador3 <- function(m) { mean(m) }
estimador4 <- function(m) { (min(m) + max(m)) / 2 }

estimarMuestra <- function(data, n_muestra, lam) {
  teta <- 1 / lam
  muestra <- data[sample(nrow(df), size=n_muestra), ]
  t1 <- apply(muestra, 1, estimador1)
  t2 <- apply(muestra, 1, estimador2)
  t3 <- apply(muestra, 1, estimador3)
  t4 <- apply(muestra, 1, estimador4)
  return(data.frame(t1, t2, t3, t4))
}
calcularMetricas <- function(datos, n_muestra, teta) {
  media <- apply(datos, 2, mean)
  varianza <- apply(datos, 2, var)
  sesgo <- media - teta
  return (data.frame(media = media, varianza = varianza, n=n_muestra, sesgo=sesgo))
}

obtenerGrafica <- function(datos, n_muestra, lam) {
  teta <- 1 / lam
  d <- pivot_longer(datos, cols = c(t1, t2, t3, t4), names_to = "categoria", values_to = "valor")
  print(d)
  bp <- ggplot(d, aes(x = categoria, y = valor, fill = categoria)) +
    geom_boxplot() +
    geom_hline(yintercept = teta, color = "red", linetype = "dashed", size = 1) +
    scale_fill_manual(values = c("#377eb8", "#ff7f00", "#4daf4a", "#f781bf")) +
    labs(title = paste("n =", n_muestra), x = "Categoría", y = "Valor") +
    theme_minimal() 
  bp <- ggplotly(bp)
  return (bp)
}
d <- estimarMuestra(df, 20, lam)
calculos <- rbind(data.frame(), calcularMetricas(d, 20, lam))
obtenerGrafica(d, 20, lam)
## # A tibble: 80 × 2
##    categoria valor
##    <chr>     <dbl>
##  1 t1        2.91 
##  2 t2        5.79 
##  3 t3        3.18 
##  4 t4        3.33 
##  5 t1        0.791
##  6 t2        1.31 
##  7 t3        1.06 
##  8 t4        1.48 
##  9 t1        2.59 
## 10 t2        4.86 
## # ℹ 70 more rows
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
d <- estimarMuestra(df, 50, lam)
calculos <- rbind(calculos, calcularMetricas(d, 50, lam))
obtenerGrafica(d, 50, lam)
## # A tibble: 200 × 2
##    categoria valor
##    <chr>     <dbl>
##  1 t1        3.25 
##  2 t2        7.59 
##  3 t3        2.72 
##  4 t4        4.07 
##  5 t1        1.47 
##  6 t2        2.18 
##  7 t3        2.03 
##  8 t4        3.18 
##  9 t1        0.418
## 10 t2        0.904
## # ℹ 190 more rows
d <- estimarMuestra(df, 100, lam)
calculos <- rbind(calculos, calcularMetricas(d, 100, lam))
obtenerGrafica(d, 100, lam)
## # A tibble: 400 × 2
##    categoria valor
##    <chr>     <dbl>
##  1 t1         4.63
##  2 t2         9.51
##  3 t3         4.59
##  4 t4         4.56
##  5 t1         1.22
##  6 t2         2.23
##  7 t3         1.35
##  8 t4         1.14
##  9 t1         2.16
## 10 t2         4.24
## # ℹ 390 more rows
d <- estimarMuestra(df, 1000, lam)
calculos <- rbind(calculos, calcularMetricas(d, 1000, lam))
obtenerGrafica(d, 1000, lam)
## # A tibble: 4,000 × 2
##    categoria valor
##    <chr>     <dbl>
##  1 t1        1.32 
##  2 t2        2.48 
##  3 t3        1.05 
##  4 t4        1.58 
##  5 t1        2.59 
##  6 t2        5.18 
##  7 t3        2.46 
##  8 t4        2.61 
##  9 t1        0.608
## 10 t2        1.13 
## # ℹ 3,990 more rows